
#####
## Responsiveness based on raw data shows little evidence of responsiveness
#####

library(texreg)
library(tidyverse)
library(reshape2)
library(estimatr)

lag.new <- function(x, n = 1L, along_with){
    index <- match(along_with - n, along_with, incomparable = NA)
    out <- x[index]
    attributes(out) <- attributes(x)
    out
}

lead.new <- function(x, n = 1L, along_with){
    index <- match(along_with + n, along_with, incomparable = NA)
    out <- x[index]
    attributes(out) <- attributes(x)
    out
}

ssm <- read.csv(file = "w_gayrights_civilunions_marriage.csv") %>%
    filter(abb != "") %>%
    mutate_at(vars(-abb), funs(as.integer(as.character(.)))) %>%
    melt(id.var = "abb", value.name = "policy") %>%
    rename(year = variable) %>%
    mutate(samesex_marriage =
               ordered(policy, labels = c("none", "civil unions", "marriage")),
           year = as.integer(str_remove(year, "^X")),
           policy_name = "same-sex marriage")

model_estimates <- read.csv("gaymarriage_disaggregated.csv")
model_estimates$abb <- gsub("abb", "", model_estimates$abb)
ssm <- merge(model_estimates, ssm,
             by.x=c("abb", "year"),
             by.y=c("abb", "year"),
             all.x=TRUE) %>%
    select(-X)

summary(lm(policy ~ naes2004_stategaymarriage_unweighted,
           data= subset(ssm, year==2012)))
summary(lm(policy ~ naes2004_stategaymarriage_unweighted + factor(abb) +
               factor(year), data= ssm))
summary(lm(policy ~ naes2004_stategaymarriage_weighted,
           data = subset(ssm, year==2012)))
reg1 <- (lm(policy ~ naes2004_stategaymarriage_weighted + factor(abb) +
                factor(year), data= ssm))


#####
## Responsiveness based on DGIRT shows evidence of responsivess
#####

require(foreign)
require(nnet)

model_estimates <- read.csv("gaymarriage_dyn_psraking8k_dgirtout.csv")
model_estimates <- rename(model_estimates, median=value_mean)
model_estimates$abb <- gsub("abb", "",model_estimates$abb)
ssm <- merge(model_estimates, ssm,  by.x=c("abb", "year"),
             by.y=c("abb", "year"), all.x = TRUE)

summary(lm(policy~median,data= subset(ssm, year==2012)))
summary(lm(policy~median + factor(abb) + factor(year),data= ssm))
ssm$policy <- as.numeric(as.vector(ssm$policy))

cor(ssm$median, ssm$naes2004_stategaymarriage_weighted, use="complete.obs")
reg2 <- (lm(policy~median + factor(abb) + factor(year),data= ssm))

summary(MASS::polr(samesex_marriage ~ median, data = subset(ssm)))

summary( multinom(policy~median + factor(abb) + factor(year), data = subset(ssm)))

stargazer(reg1, reg2)


### table 2
fit_responsiveness_model <- function (data, xvar, state_fe = TRUE) {
    if (state_fe)
        f <- paste("y ~ ", xvar, "+ abb + Year")
    else
        f <- paste("y ~ ", xvar, "+ Year")
    print(f)
    data[xvar] <- data[xvar] * 100
    mod <- data %>%
        mutate(Year = factor(year),
               y = as.numeric(samesex_marriage)) %>% 
        lm_robust(as.formula(f), data = ., clusters = abb)
    mod
}

ssm %>%
    fit_responsiveness_model("median") 
 

xvars <- c("naes2004_stategaymarriage_weighted", "median")

summary(ssm[xvars])
cor(ssm[xvars], use = "complete.obs")

reg_dyn.ls <- plyr::llply(xvars, function (xv) {
    ssm %>%
        fit_responsiveness_model(data = ., xvar = xv, state_fe = TRUE)
})

reg_xs.ls <- plyr::llply(xvars, function (xv) {
    ssm %>%
        fit_responsiveness_model(data = ., xvar = xv, state_fe = FALSE)
})

custom.ls <- as.list(rep("Mass Support for Same-Sex Marriage (%)", 2))
names(custom.ls) <- xvars

texreg(c(reg_xs.ls, reg_dyn.ls), dcolumn = TRUE, include.ci = FALSE,
       custom.coef.map = custom.ls,
       digits = 3, use.packages = FALSE, include.adjrs = FALSE,
       include.rmse = FALSE, table = FALSE,
       custom.model.names = c("Disaggregated (CS)",
                              "Smoothed (CS)",
                              "Disaggregated (TSCS)",
                              "Smoothed (TSCS)"))

?texreg
